﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// Provides a method to deal with polynomial interpolations by the theorem of Neville.
    /// </summary>
    [Serializable]
    public static class NevillePolynomialInterpolation
    {
        /// <summary>
        /// Polynomialinterpolation by the theorem of Neville.
        /// </summary>
        /// <param name="interpolatingPoints">The number of interpolating points.</param>
        /// <param name="x">The array, which containing interpolating points for x.</param>
        /// <param name="y">The array, which containing interpolating points for f(x).</param>
        /// <param name="t">The value for evaluation.</param>
        /// <returns>The value f(x) at the position x.</returns>
        public static double NevilleInterpolation(int interpolatingPoints, double[] x, double[] y, double t)
        {
            for (int i = 1; i < interpolatingPoints; i++)
            {
                for (int j = interpolatingPoints - 1; j >= i; j--)
                {
                    y[j] = ((t - x[j - i]) * y[j] - (t - x[j]) * y[j - 1]) / (x[j] - x[j - i]);
                }
            }

            return y[interpolatingPoints - 1];
        }

        /// <summary>
        /// Polynomialinterpolation by the theorem of Neville at position x=0.
        /// </summary>
        /// <param name="interpolatingPoints">The number of interpolating points.</param>
        /// <param name="x">The array, which containing interpolating points for x.</param>
        /// <param name="y">The array, which containing interpolating points for f(x).</param>
        /// <returns>The value f(x) at the position x.</returns>
        public static double NevilleInterpolation(int interpolatingPoints, double[] x, double[] y)
        {
            for (int i = 1; i < interpolatingPoints; i++)
            {
                for (int j = interpolatingPoints - 1; j >= i; j--)
                {
                    y[j] = ((- x[j - i]) * y[j] - (- x[j]) * y[j - 1]) / (x[j] - x[j - i]);
                }
            }

            return y[interpolatingPoints - 1];
        }
    }
}